## Replication code for 
## "An Experimental Test of the Effects of Fear in a Coordination Game"
## Journal of Experimental Political Science 
## Aldama, Sambrano, Vásquez-Cortés, Young

#####
## Prepare and clean data
#####
 
rm(list = ls())
library(zTree)
library(tidyverse)
library(data.table)
library(stargazer)
library(mediation)
library(lmtest)
library(sandwich)
library("Hmisc")

setwd("~/Dropbox/Fear and coordination/11_JEPS_Replication_Files/")


#####
## Cleaning + building dataset
#####

# Read in raw data by laboratory session
zTT <- zTreeTables(c("01_Data/New data/190418_1853.xls", "01_Data/New data/190425_1928.xls",
                     "01_Data/New data/190502_1923.xls", "01_Data/New data/190510_1658.xls",
                     "01_Data/New data/190514_1708.xls", "01_Data/New data/190515_1703.xls",
                     "01_Data/New data/190607_1627.xls", "01_Data/New data/190612_1631.xls",
                     "01_Data/New data/190613_1617.xls", "01_Data/New data/190614_1611.xls",
                     "01_Data/New data/191122_1119 (version 3).xls",
                     "01_Data/New data/191122_1013.xls", 
                     "01_Data/New data/191122_1238.xls", "01_Data/New data/191122_1435.xls",
                     "01_Data/New data/191125_0923.xls", "01_Data/New data/191125_1101.xls",
                     "01_Data/New data/191125_1357.xls", "01_Data/New data/191126_0833.xls",
                     "01_Data/New data/191126_1058.xls", "01_Data/New data/191126_1305.xls",
                     "01_Data/New data/191126_1438.xls", "01_Data/New data/191202_1014.xls",
                     "01_Data/New data/191202_1112.xls", "01_Data/New data/191202_1225.xls",
                     "01_Data/New data/191202_1355.xls", "01_Data/New data/191202_1534.xls",
                     "01_Data/New data/191203_1053.xls", "01_Data/New data/191203_1228.xls",
                     "01_Data/New data/191203_1356.xls", "01_Data/New data/191206_0906.xls",
                     "01_Data/New data/191206_1325.xls", "01_Data/New data/191206_1400.xls"))

with(zTT$subjects, table(Treatment,Period))
dat <- data.frame(zTT$subjects)
dat$subject.id <- paste0(dat$Subject, 'day', dat$Date)
dat$Set <- ifelse(dat$Period <= 5, 1,
                  ifelse(dat$Period > 5 & dat$Period <= 10, 2,
                         ifelse(dat$Period > 10 & dat$Period <= 15, 3, 4)))
dat$Treatment <- NULL

## globals (Y and C)
glb <- zTT$globals
glb <- subset(glb, select = c(Date, Period, y, c))
dat <- merge(dat, glb, by = c("Date", "Period"), all.x = T)

## Data with questionnaire (covariates)
sbj <- zTreeSbj(c('01_Data/New data/190418_1853.sbj', '01_Data/New data/190425_1928.sbj',
                  '01_Data/New data/190502_1923.sbj', '01_Data/New data/190510_1658.sbj',
                  '01_Data/New data/190514_1708.sbj', '01_Data/New data/190515_1703.sbj',
                  '01_Data/New data/190607_1627.sbj', '01_Data/New data/190612_1631.sbj',
                  '01_Data/New data/190613_1617.sbj', '01_Data/New data/190614_1611.sbj', 
                  "01_Data/New data/191122_1119.sbj", "01_Data/New data/191122_1013.sbj", 
                  "01_Data/New data/191122_1238.sbj", "01_Data/New data/191122_1435.sbj",
                  "01_Data/New data/191125_0923.sbj", "01_Data/New data/191125_1101.sbj",
                  "01_Data/New data/191125_1357.sbj", "01_Data/New data/191126_0833.sbj",
                  "01_Data/New data/191126_1058.sbj", "01_Data/New data/191126_1305.sbj",
                  "01_Data/New data/191126_1438.sbj", "01_Data/New data/191202_1014.sbj",
                  "01_Data/New data/191202_1112.sbj", "01_Data/New data/191202_1225.sbj",
                  "01_Data/New data/191202_1355.sbj", "01_Data/New data/191202_1534.sbj",
                  "01_Data/New data/191203_1053.sbj", "01_Data/New data/191203_1228.sbj",
                  "01_Data/New data/191203_1356.sbj", "01_Data/New data/191206_0906.sbj",
                  "01_Data/New data/191206_1325.sbj", "01_Data/New data/191206_1400.sbj"))

sbj$subject.id <- paste0(sbj$Subject, 'day', sbj$Date)
sbj$female <- ifelse(sbj$sex=="1 = Male;", 0, 1)
sbj$age <- as.numeric(as.character(sbj$age))
sbj$signal_bias <- as.numeric(as.character(strtrim(sbj$perceptionsignal, 1)))-2
sbj$partner_cautious <- as.numeric(as.character(strtrim(sbj$cautious, 1)))-2
sbj$owncutoff <- as.numeric(as.character(sbj$owncutoff))
sbj$othercutoff <- as.numeric(as.character(sbj$othercutoff))
sbj$Date <- NULL; sbj$Subject <- NULL
sbj$fear <- as.numeric(strtrim(sbj$fear, 1))
sbj$anger <- as.numeric(strtrim(sbj$anger, 1))
sbj$disgust <- as.numeric(strtrim(sbj$disgust, 1))
sbj$sadness <- as.numeric(strtrim(sbj$sadness, 1))
sbj$happiness <- as.numeric(strtrim(sbj$happiness, 1))
sbj$surprise <- as.numeric(strtrim(sbj$surprise, 1))
sbj$computer <- gsub("Lab|lab", "", sbj$client)

sbj$baselinetube2 <- as.numeric(strtrim(sbj$baselinetube, 3))
sbj$treatmenttesttube2 <- as.numeric(strtrim(sbj$treatmenttesttube, 3))

# ## Export clean data
# write.dta(sbj, "~/sbj.dta")

dat <- merge(dat, sbj, by = 'subject.id', all.x = T)

## Data with randomization info: indicate types of rounds
rand <- read.csv('01_Data/New data/randomization_new.csv')
dat <- merge(dat, rand, by = c('Date', 'Set'), all.x = T)

## collapse to each subject 
dat2 <- data.table(dat)
dat2 <- dat2[, list('coop_pre' = mean(action_i[Type=="pre"]),
                    'coop_avg' = mean(action_i[Type!="pre"]),
                    'coop_h1n1' = mean(action_i[Type=="h1n1"]),
                    'coop_h1n0' = mean(action_i[Type=="h1n0"]),
                    'coop_h0n1' = mean(action_i[Type=="h0n1"]),
                    'signal_h1n1' = mean(signal[Type=='h1n1']),
                    'signal_h1n0' = mean(signal[Type=='h1n0']),
                    'signal_h0n1' = mean(signal[Type=='h0n1']),
                    'order_h1n1' = mean(Set[Type=='h1n1'])-1,
                    'order_h1n0' = mean(Set[Type=='h1n0'])-1,
                    'order_h0n1' = mean(Set[Type=='h0n1'])-1),
             by = c('subject.id')]

## # Export Data to analysis in Stata
#write.csv(dat2, '~/data_clean_new_notcollapse.csv')

## calculate spread of round
dat$spread <- dat$y + dat$c
table(dat$spread, dat$Treatment)
table(dat$spread, dat$Treatment, dat$Type)
table(dat$spread, dat$Type)
dat$spread_st <- (dat$spread - mean(dat$spread, na.rm=T))/sd(dat$spread, na.rm=T)

## calculate pre-treatment average cooperation
dat <- data.table(dat)
dat[, 'pre_avg' := mean(action_i[Type=="pre"]), by = 'subject.id']

#####
## Data for analysis
#####

dat3 <- join(dat, dat2, by='subject.id', type='left', match='all')
#save(dat3,file = "~/dat3.RData")

#####
## FIGURE 3: Effect of Fear on Mobilization
#####

#load("dat3.RData")

data_h1n1 <- subset(dat3, Type=="h1n1")
mod1a <- lm(action_i ~ Treatment, data = data_h1n1)
mod1a$se <- coeftest(mod1a, vcov = vcovCL, cluster = ~Date)[,2]
mod1b <- lm(action_i ~ Treatment + female + signal + age + pre_avg, data = data_h1n1)
mod1b$se <- coeftest(mod1b, vcov = vcovCL, cluster = ~Date)[,2]

data_h1n0 <- subset(dat3, Type=="h1n0")
mod2a <- lm(action_i ~ Treatment, data = data_h1n0)
mod2a$se <- coeftest(mod2a, vcov = vcovCL, cluster = ~Date)[,2]
mod2b <- lm(action_i ~ Treatment + female + signal + age + pre_avg, data = data_h1n0)
mod2b$se <- coeftest(mod2b, vcov = vcovCL, cluster = ~Date)[,2]

data_h0n1 <- subset(dat3, Type=="h0n1")
mod3a <- lm(action_i ~ Treatment, data = data_h0n1)
mod3a$se <- coeftest(mod3a, vcov = vcovCL, cluster = ~Date)[,2]
mod3b <- lm(action_i ~ Treatment + female + signal + age + pre_avg, data = data_h0n1)
mod3b$se <- coeftest(mod3b, vcov = vcovCL, cluster = ~Date)[,2]

coefs <- c(summary(mod1a)$coef['Treatment',1], 
           summary(mod2a)$coef['Treatment',1], 
           summary(mod3a)$coef['Treatment',1])
ses <- c(mod1a$se['Treatment'], 
         mod2a$se['Treatment'],
         mod3a$se['Treatment'])

plotdat <- data.frame('Treatment' = rep(c('Treatment'),3),
                      x = c('Standard', 'Full Information', 'Computer'),
                      coefs, ses, 
                      'up95' = coefs + 1.95*ses, 'lo95' = coefs - 1.95*ses)

ggplot(plotdat) +
  geom_hline(yintercept = 0, colour = gray(1/2), lty = 2) +
  geom_errorbar(aes(x = x, y = coefs, ymin = lo95, ymax = up95),
                lwd = 1, position = position_dodge(width = 1/2),
                width = 0.25) +
  geom_point(aes(x = x, y = coefs, ymin = lo95, ymax = up95), 
             position = position_dodge(width=1/2)) +
  theme_minimal() +
  theme(text = element_text(size = 20),
        axis.text.x = element_text(angle = 45, hjust = 1)) +
  scale_x_discrete(name ="Type of Round", 
                   limits=c("Standard","Full Information", 'Computer')) +
  scale_y_continuous(name = 'Treatment Effect')




#####
## TABLE 3: Differential Effect of Fear by Types of Rounds
#####

dat3 <- dat3 %>% 
  mutate('computer' = if_else(Type=="h0n1", 1, 0),
         'noiseless' = if_else(Type=="h1n0", 1, 0))

dat3_compX <- dat3 %>% filter(Type == "h0n1" | Type == "h1n1")
dat3_noiselessX <- dat3 %>% filter(Type == "h1n0" | Type == "h1n1")

mod1a <- lm(action_i ~ Treatment + computer + Treatment:computer, data = dat3_compX)
mod1a$se <- coeftest(mod1a, vcov = vcovCL, cluster = ~Date)[,2]
mod1b <- lm(action_i ~ Treatment + computer + Treatment:computer + signal, data = dat3_compX)
mod1b$se <- coeftest(mod1b, vcov = vcovCL, cluster = ~Date)[,2]
mod1c <- lm(action_i ~ Treatment + computer + Treatment:computer + female + age + signal + pre_avg, data = dat3_compX)
mod1c$se <- coeftest(mod1c, vcov = vcovCL, cluster = ~Date)[,2]

mod2a <- lm(action_i ~ Treatment + noiseless + Treatment:noiseless, data = dat3_noiselessX)
mod2a$se <- coeftest(mod2a, vcov = vcovCL, cluster = ~Date)[,2]
mod2b <- lm(action_i ~ Treatment + noiseless + Treatment:noiseless + signal, data = dat3_noiselessX)
mod2b$se <- coeftest(mod2b, vcov = vcovCL, cluster = ~Date)[,2]
mod2c <- lm(action_i ~ Treatment + noiseless + Treatment:noiseless + female + age + signal + pre_avg, data = dat3_noiselessX)
mod2c$se <- coeftest(mod2c, vcov = vcovCL, cluster = ~Date)[,2]

stargazer(mod1a, mod1b, mod1c, mod2a, mod2b, mod2c,
          se = list(mod1a$se, mod1b$se, mod1c$se, mod2a$se, mod2b$se, mod2c$se),
          no.space = T, 
          covariate.labels = c('Fear Treatment', 
                               'Computer Round', 
                               'Female', 'Age', 'Signal',
                               'Pre-Treatment Cooperation',
                               'Treatment \\times Computer Round',
                               'Noiseless Round',
                               'Treatment \\times Noiseless Round',
                               'Intercept'),
          keep.stat = c('n', 'rsq'))


#####
## FIGURE 5: Effect of Treatment on Emotions
#####
## emotions
mod1 <- lm(fear ~ Treatment, data = dat3)
mod1$se <- coeftest(mod1, vcov = vcovCL, cluster = ~Date)[,2]

mod2 <- lm(anger ~ Treatment , data = dat3)
mod2$se <- coeftest(mod2, vcov = vcovCL, cluster = ~Date)[,2]

mod3 <- lm(sadness ~ Treatment , data = dat3)
mod3$se <- coeftest(mod3, vcov = vcovCL, cluster = ~Date)[,2]

mod4 <- lm(surprise ~ Treatment , data = dat3)
mod4$se <- coeftest(mod4, vcov = vcovCL, cluster = ~Date)[,2]

mod5 <- lm(happiness ~ Treatment , data = dat3)
mod5$se <- coeftest(mod5, vcov = vcovCL, cluster = ~Date)[,2]

mod6 <- lm(disgust ~ Treatment , data = dat3)
mod6$se <- coeftest(mod6, vcov = vcovCL, cluster = ~Date)[,2]

coefs <- c(summary(mod1)$coef['Treatment',1], 
           summary(mod2)$coef['Treatment',1], 
           summary(mod3)$coef['Treatment',1], 
           summary(mod4)$coef['Treatment',1], 
           summary(mod5)$coef['Treatment',1], 
           summary(mod6)$coef['Treatment',1])
ses <- c(mod1$se['Treatment'], 
         mod2$se['Treatment'],
         mod3$se['Treatment'],
         mod4$se['Treatment'],
         mod5$se['Treatment'],
         mod6$se['Treatment'])

plotdat <- data.frame('Treatment' = rep(c('Treatment'),6),
                      x = c('Fear', 'Anger', 'Sadness', 'Surprise', 
                            'Happiness', 'Disgust'),
                      coefs, ses, 
                      'up95' = coefs + 1.96*ses, 'lo95' = coefs - 1.96*ses)

ggplot(plotdat) +
  geom_hline(yintercept = 0, colour = gray(1/2), lty = 2) +
  geom_errorbar(aes(x = x, y = coefs, ymin = lo95, ymax = up95),
                lwd = 1, position = position_dodge(width = 1/2),
                width = 0.25) +
  geom_point(aes(x = x, y = coefs, ymin = lo95, ymax = up95), 
             position = position_dodge(width=1/2)) +
  theme_minimal() +
  theme(text = element_text(size = 20),
        axis.text.x = element_text(angle = 45, hjust = 1)) +
  scale_x_discrete(name ="Emotion") +
  scale_y_continuous(name = 'Treatment Effect')



#####
## TABLE 5 (APPENDIX): Effect of fear treatment in Standard rounds
#####

# data_h1n1
data_h1n1 <- subset(dat3, Type=="h1n1")
mod1a <- lm(action_i ~ Treatment, data = data_h1n1)
mod1a$se <- coeftest(mod1a, vcov = vcovCL, cluster = ~Date)[,2]

mod1b <- lm(action_i ~ Treatment + order_h1n1 + signal, data = data_h1n1)
mod1b$se <- coeftest(mod1b, vcov = vcovCL, cluster = ~Date)[,2]

mod1c <- lm(action_i ~ Treatment + order_h1n1 + signal + coop_pre + female + age, data = data_h1n1)
mod1c$se <- coeftest(mod1c, vcov = vcovCL, cluster = ~Date)[,2]

stargazer(mod1a, mod1b, mod1c,
          se = list(mod1a$se, mod1b$se, mod1c$se),
          no.space = T, 
          covariate.labels = c('Fear Treatment', 'Order of Round',
                               'Signal', 'Pre-Treatment Cooperation',
                               'Female', 'Age', 'Intercept'),
          keep.stat = c('n', 'rsq'))

#####
## TABLE 6 (APPENDIX): Effect of fear treatment in Full Information Rounds
#####
# data_h1n0
mod2a <- lm(action_i ~ Treatment, data = data_h1n0)
mod2a$se <- coeftest(mod2a, vcov = vcovCL, cluster = ~Date)[,2]

mod2b <- lm(action_i ~ Treatment + order_h1n0 + signal, data = data_h1n0)
mod2b$se <- coeftest(mod2b, vcov = vcovCL, cluster = ~Date)[,2]

mod2c <- lm(action_i ~ Treatment + order_h1n0 + signal + coop_pre + female + age, data = data_h1n0)
mod2c$se <- coeftest(mod2c, vcov = vcovCL, cluster = ~Date)[,2]

stargazer(mod2a, mod2b, mod2c,
          se = list(mod2a$se, mod2b$se, mod2c$se),
          no.space = T, 
          covariate.labels = c('Fear Treatment', 'Order of Round',
                               'Signal', 'Pre-Treatment Cooperation',
                               'Female', 'Age', 'Intercept'),
          keep.stat = c('n', 'rsq'))


#####
## TABLE 7 (APPENDIX): Effect of fear treatment in Computer Rounds
#####

# data_h0n1
mod3a <- lm(action_i ~ Treatment, data = data_h0n1)
mod3a$se <- coeftest(mod3a, vcov = vcovCL, cluster = ~Date)[,2]

mod3b <- lm(action_i ~ Treatment + order_h0n1 + signal, data = data_h0n1)
mod3b$se <- coeftest(mod3b, vcov = vcovCL, cluster = ~Date)[,2]

mod3c <- lm(action_i ~ Treatment + order_h0n1 + signal + coop_pre + female + age, data = data_h0n1)
mod3c$se <- coeftest(mod3c, vcov = vcovCL, cluster = ~Date)[,2]

stargazer(mod3a, mod3b, mod3c,
          se = list(mod3a$se, mod3b$se, mod3c$se),
          no.space = T, 
          covariate.labels = c('Fear Treatment', 'Order of Round',
                               'Signal', 'Pre-Treatment Cooperation',
                               'Female', 'Age', 'Intercept'),
          keep.stat = c('n', 'rsq'))


#####
## TABLE 8 (APPENDIX): Effect of fear treatment by payoff spread
#####

d <- dat[dat3$Type=="h1n1",]
mod4a <- lm(action_i ~ Treatment + spread + Treatment:spread, data = d)
mod4a$se <- coeftest(mod4a, vcov = vcovCL, cluster = ~Date)[,2]

d2 <- dat[dat3$Type=="h0n1",]
mod4a2 <- lm(action_i ~ Treatment + spread + Treatment:spread, data = d2)
mod4a2$se <- coeftest(mod4a2, vcov = vcovCL, cluster = ~Date)[,2]

d3 <- dat[dat3$Type=="h1n0",]
mod4a3 <- lm(action_i ~ Treatment + spread + Treatment:spread, data = d3)
mod4a3$se <- coeftest(mod4a3, vcov = vcovCL, cluster = ~Date)[,2]

stargazer(mod4a, mod4a2, mod4a3, 
          se = list(mod4a$se, mod4a2$se, mod4a3$se),
          no.space = T,
          keep.stat = c('n', 'rsq'))



#####
## Figure 6 (APPENDIX):  Probability of Mobilization by Round
#####

data_h1n1$Fear <- factor(data_h1n1$Treatment)
data_h1n0$Fear <- factor(data_h1n0$Treatment)
data_h0n1$Fear <- factor(data_h0n1$Treatment)

aa <- ggplot(data_h1n1, aes(x=signal, y=action_i, colour=Fear)) +
  geom_point(position = position_jitter(width = 0.3, height = 0.06), alpha=0.4, shape=21, size=1.5)+
  stat_smooth(method = "glm", method.args = list(family="binomial"), se=FALSE)+
  geom_vline(xintercept = 1, colour = gray(1/2), lty = 2)+
  ggtitle("Standard")+
  scale_colour_discrete(labels = c("Control", "Treatment")) +
  theme(text = element_text(size = 15))+
  scale_y_continuous(name = "Probability of Mobilization", breaks = c(0,0.5,1)) + 
  scale_x_continuous(name = "Signal")

bb <- ggplot(data_h1n0, aes(x=signal, y=action_i, colour=Fear)) +
  geom_point(position = position_jitter(width = 0.3, height = 0.06), alpha=0.4, shape=21, size=1.5)+
  stat_smooth(method = "glm", method.args = list(family="binomial"), se=FALSE)+
  geom_vline(xintercept = 1, colour = gray(1/2), lty = 2)+
  ggtitle("Full Information")+
  scale_colour_discrete(labels = c("Control", "Treatment")) +
  theme(text = element_text(size = 15))+
  scale_y_continuous(name = "Probability of Mobilization", breaks = c(0,0.5,1)) +  
  scale_x_continuous(name = "Signal")

cc <- ggplot(data_h0n1, aes(x=signal, y=action_i, colour=Fear)) +
  geom_point(position = position_jitter(width = 0.3, height = 0.06), alpha=0.4, shape=21, size=1.5) +
  stat_smooth(method = "glm", method.args = list(family="binomial"), se=FALSE) +
  geom_vline(xintercept = 1, colour = gray(1/2), lty = 2)+
  ggtitle("Computer") +
  scale_colour_discrete(labels = c("Control", "Treatment")) +
  theme(text = element_text(size = 15))+
  scale_y_continuous(name = "Probability of Mobilization", breaks = c(0,0.5,1)) +  
  scale_x_continuous(name = "Signal")

combined <- aa + bb + cc  & theme(legend.position = "bottom")
combined + plot_layout(guides = "collect")


#####
## Table 9 (APPENDIX):  Estimated Threshold by Type of Round
#####

data_h1n1T <- subset(data_h1n1, Treatment=="1")
data_h1n1C <- subset(data_h1n1, Treatment=="0")

data_h1n0T <- subset(data_h1n0, Treatment=="1")
data_h1n0C <- subset(data_h1n0, Treatment=="0")

data_h0n1T <- subset(data_h0n1, Treatment=="1")
data_h0n1C <- subset(data_h0n1, Treatment=="0")

logith1n1T <- glm(action_i ~ signal, data = data_h1n1T, family = "binomial")
Thresholdh1n1T <- -coef(logith1n1T)["(Intercept)"]/summary(logith1n1T)$coef['signal',1]
sdThresh1n1T <- -pi/summary(logith1n1T)$coef['signal',1]*sqrt(3)

logith1n1C <- glm(action_i ~ signal, data = data_h1n1C, family = "binomial")
Thresholdh1n1C <- -coef(logith1n1C)["(Intercept)"]/summary(logith1n1C)$coef['signal',1]
sdThresh1n1C <- -pi/summary(logith1n1C)$coef['signal',1]*sqrt(3)

logith1n0T <- glm(action_i ~ signal , data = data_h1n0T, family = "binomial")
Thresholdh1n0T <- - coef(logith1n0T)["(Intercept)"]/summary(logith1n0T)$coef['signal',1]
sdThresh1n0T <- -pi/summary(logith1n0T)$coef['signal',1]*sqrt(3)

logith1n0C <- glm(action_i ~ signal , data = data_h1n0C, family = "binomial")
Thresholdh1n0C <- - coef(logith1n0C)["(Intercept)"]/summary(logith1n0C)$coef['signal',1]
sdThresh1n0C <- -pi/summary(logith1n0C)$coef['signal',1]*sqrt(3)

logith0n1T <- glm(action_i ~ signal , data = data_h0n1T, family = "binomial")
Thresholdh0n1T <- - coef(logith0n1T)["(Intercept)"]/summary(logith0n1T)$coef['signal',1]
sdThresh0n1T <- -pi/summary(logith0n1T)$coef['signal',1]*sqrt(3)

logith0n1C <- glm(action_i ~ signal , data = data_h0n1C, family = "binomial")
Thresholdh0n1C <- - coef(logith0n1C)["(Intercept)"]/summary(logith0n1C)$coef['signal',1]
sdThresh0n1C <- -pi/summary(logith0n1C)$coef['signal',1]*sqrt(3)

threshold_table <- rbind(c(Thresholdh1n1T, Thresholdh0n1T),
                         c(Thresholdh1n1C, Thresholdh0n1C))
colnames(threshold_table) <- rep(c('Standard', 'Computer'))
rownames(threshold_table) <- c('Fear', 'Control')
xtable(threshold_table)


#####
## Figure 7 (APPENDIX):  Rational choices by type of round - STATA
#####

## In stata do file:  Figure7.do

#####
## Table 10 (APPENDIX): Effect of Treatment on levels of Alpha-Amylase
#####

amylase <- read_excel("01_Data/New data/survey with alpha amylase.xlsx")

dat4 <- merge(dat3, amylase, by.x = 'subject.id', by.y = 'subject_id')
dat4 <- dat4[ which(dat4$numberappearsmultipletimes==0),]

## collapse by subject
dat5 <- dat4[, list('amylase_base' = mean(amylasebaseline),
                    'amylase_treatment' = mean(amylasetreatment),
                    'Treatment' = mean(Treatment)),
             by = c('subject.id', 'Date')]

amy1 <- lm(amylase_treatment ~ Treatment + amylase_base, data = dat5)
amy1$se <- coeftest(amy1, vcov = vcovCL, cluster = ~Date)[,2]

stargazer(amy1,
          se = list(amy1$se),
          no.space = T,
          covariate.labels = c('Fear Treatment', 'Amylase baseline'),
          keep.stat = c('n', 'rsq'))

#####
## Table 11 (APPENDIX): Alpha-Amylase as Predictor
#####

amylase <- read_excel("01_Data/New data/survey with alpha amylase.xlsx")

dat4 <- merge(dat3, amylase, by.x = 'subject.id', by.y = 'subject_id')

## collapse by subject
dat5 <- dat4[, list('amylase_base' = mean(amylasebaseline),
                    'amylase_treatment' = mean(amylasetreatment),
                    'Treatment' = mean(Treatment),
                    'Mobilization' = mean(action_i),
                    'signal' = mean(signal)),
             by = c('subject.id', 'Date')]

dat6 <- merge(dat5, sbj, by = 'subject.id', all.x = T)
dat6$amylase_diff <- dat6$amylase_treatment - dat6$amylase_base

mod_Amy <- lm(Mobilization ~ amylase_diff, data = dat6)
mod_Amy$se <- coeftest(mod_Amy, vcov = vcovCL, cluster = ~Date)[,2]

mod_Amy2 <- lm(Mobilization ~ amylase_diff + signal, data = dat6)
mod_Amy2$se <- coeftest(mod_Amy2, vcov = vcovCL, cluster = ~Date)[,2]

mod_Amy3 <- lm(Mobilization ~ amylase_diff + signal + female + age, data = dat6)
mod_Amy3$se <- coeftest(mod_Amy3, vcov = vcovCL, cluster = ~Date)[,2]

stargazer(mod_Amy, mod_Amy2, mod_Amy3,
          se = list(mod_Amy$se, mod_Amy2$se, mod_Amy3$se),
          no.space = T, 
          covariate.labels = c('Amylase Difference',
                               'Signal',
                               'Female', 'Age', 'Intercept'),
          keep.stat = c('n', 'rsq'))



#####
## FIGURE 10: Effect of Fear on Mobilization - First Five Rounds
#####

data_F5 <- subset(dat3, Period>5 & Period<11)

data_h1n1_F5 <- subset(data_F5, Type=="h1n1")
mod1a_F5 <- lm(action_i ~ Treatment, data = data_h1n1_F5)
mod1a_F5$se <- coeftest(mod1a_F5, vcov = vcovCL, cluster = ~Date)[,2]

data_h1n0_F5 <- subset(data_F5, Type=="h1n0")
mod2a_F5 <- lm(action_i ~ Treatment, data = data_h1n0_F5)
mod2a_F5$se <- coeftest(mod2a_F5, vcov = vcovCL, cluster = ~Date)[,2]

data_h0n1_F5 <- subset(data_F5, Type=="h0n1")
mod3a_F5 <- lm(action_i ~ Treatment, data = data_h0n1_F5)
mod3a_F5$se <- coeftest(mod3a_F5, vcov = vcovCL, cluster = ~Date)[,2]

coefs <- c(summary(mod1a_F5)$coef['Treatment',1], 
           summary(mod2a_F5)$coef['Treatment',1], 
           summary(mod3a_F5)$coef['Treatment',1])
ses <- c(mod1a_F5$se['Treatment'], 
         mod2a_F5$se['Treatment'],
         mod3a_F5$se['Treatment'])

plotdat <- data.frame('Treatment' = rep(c('Treatment'),3),
                      x = c('Standard', 'Full Information', 'Computer'),
                      coefs, ses, 
                      'up95' = coefs + 1.95*ses, 'lo95' = coefs - 1.95*ses)

ggplot(plotdat) +
  geom_hline(yintercept = 0, colour = gray(1/2), lty = 2) +
  geom_errorbar(aes(x = x, y = coefs, ymin = lo95, ymax = up95),
                lwd = 1, position = position_dodge(width = 1/2),
                width = 0.25) +
  geom_point(aes(x = x, y = coefs, ymin = lo95, ymax = up95), 
             position = position_dodge(width=1/2)) +
  theme_minimal() +
  theme(text = element_text(size = 20),
        axis.text.x = element_text(angle = 45, hjust = 1)) +
  scale_x_discrete(name ="Type of Round", 
                   limits=c("Standard","Full Information", 'Computer')) +
  scale_y_continuous(name = 'Treatment Effect')


#####
## FIGURES 11, 12 and 13: Mediation analysis 
#####

# Mediation fear
subdata <- subset(dat3, Type=="h1n1")
subdata <- subset(subdata, action_i!="NA")

med.fit <- lm(fear ~ Treatment, data = subdata)
out.fit <- lm(action_i ~ Treatment + fear, data = subdata)

med.out <- mediate(med.fit, out.fit, treat = "Treatment", mediator = "fear", robustSE = TRUE, sims = 1000)
plot(med.out)

# Mediation surprise
subdata <- subset(dat3, Type=="h1n1")
subdata <- subset(subdata, action_i!="NA")

med.fit2 <- lm(surprise ~ Treatment, data = subdata)
out.fit2 <- lm(action_i ~ Treatment + surprise, data = subdata)

med.out2 <- mediate(med.fit2, out.fit2, treat = "Treatment", mediator = "surprise", robustSE = TRUE, sims = 1000)
plot(med.out2)


# Mediation both 
subdata$emotions_both <-  subdata$fear + subdata$surprise

med.fit3 <- lm(emotions_both ~ Treatment, data = subdata)
out.fit3 <- lm(action_i ~ Treatment + emotions_both, data = subdata)

med.out3 <- mediate(med.fit3, out.fit3, treat = "Treatment", mediator = "emotions_both", robustSE = TRUE, sims = 1000)
plot(med.out3)





